xxxxxxxxxx

2.7 ms

xxxxxxxxxx

83.9 ms

xxxxxxxxxx

565 s

Statistical inference with Julia

Group Retreat Heidelberg, Matthias and Maurice, 15. September 2020

xxxxxxxxxx

395 ms

Poisson parameter:

xxxxxxxxxx

9.8 ms

xxxxxxxxxx

254 ms

(Package imports above title)

xxxxxxxxxx

24.3 μs

Before we start

Tutorial requirements (the complete process should take <30 min.)

Julia installation
Optional: Atom and Juno (the Julia editor)
Julia packages listed below:

xxxxxxxxxx

1.6 ms

Status `/private/var/folders/29/l75f_tqn5fzb5qylqt0w2xbr00012v/T/jl_q64l4i/Project.toml`
  [6e4b80f9] BenchmarkTools v0.5.0
  [0c46a032] DifferentialEquations v6.15.0
  [31c24e10] Distributions v0.23.12
  [4138dd39] JLD v0.10.0
  [91a5bcdd] Plots v1.6.6
  [c3e4b0f8] Pluto v0.11.14
  [7f904dfe] PlutoUI v0.6.4
  [f3b207a7] StatsPlots v0.14.13
  [fce5fe82] Turing v0.14.4

xxxxxxxxxx

367 ms

How to install and precompile Julia packages

Start Julia in the terminal (Applications -> Julia-1.5)
Enter “]” to get into Julia package mode (pkg>)
Install a package X by writing “add X” (e.g., "add Plots", repeat for each package)
View all installed packages with “status”
Precompile all packages with “precompile”
Leave the package mode by pressing backspace (julia>)
(General info: Import a package by “using X” for a package X)

How to start this Pluto notebook

Start Julia in the terminal (Applications -> Julia-1.5)
Import Pluto by typing "using Pluto"
Open Pluto in a browser by typing "Pluto.run(1234)"
Open this notebook file via the Pluto user interface

(Currently, Pluto is supposed to work best with Firefox or Chrome)

xxxxxxxxxx

10.9 μs

Content of this Tutorial

Introduction to Julia
Performance tips and a Gillespie algorithm
Differential Equations in Julia
Bayesian Inference and Probabilistic Programming
Everything from above, just combined :)
Further Information/References

xxxxxxxxxx

5.3 μs

1. Introduction: What is Julia?

Goal: Write code as in Matlab/Python/R which is as fast as C/C++/Fortran

Basic language principles

Multi-platform via LLVM
Multi-paradigm (functional, object-oriented, ...)
Multiple dispatch
Dynamic with type inference
Just-in-time (JIT) compilation
Broad metaprogramming features
Open source

xxxxxxxxxx

10 μs

The Two Language Problem

Effort and language barriers

prototype code with a language that is fast to code with (Python, R, Matlab)
tranfer the prototype code to a language that runs fast (C, C++, Fortran)

The numpy package in Python (basic array class for fast computing)

A typical Julia package (Turing.jl)

Julia solves the two language problem

you can prototype code in Julia
you can improve code performance in the same language
no language barrier

xxxxxxxxxx

25.7 ms

Basic syntax

Getting help

Precede unknown object/expression with a question mark (?) to see the docs (Pluto has Live docs for this)

xxxxxxxxxx

5.9 μs

Some basic data structures

Tuples (ordered, immutable)
Named Tuples (ordered, immutable)
Arrays (ordered, mutable)
Dictonaries (unordered, mutable)

Tuples and Named Tuples are constructed using paratheses

xxxxxxxxxx

9.5 μs

xxxxxxxxxx
 
myTuple = (3, 1, 4);

3.1 μs

xxxxxxxxxx
 
myNamedTuple = (foo = 1, bar = 5, baz = 9);

4 ms

Elements are accessed by position (In Julia indexing starts at 1!)

xxxxxxxxxx

6.6 μs

xxxxxxxxxx
 
myTuple[2];

3.9 μs

xxxxxxxxxx
 
myNamedTuple[3];

2.8 ms

Named Tuples can additionally be accessed by name

xxxxxxxxxx

7.3 μs

xxxxxxxxxx
 
myNamedTuple[:baz];

3.1 ms

To form collections of related values, use Dictionaries or Arrays

xxxxxxxxxx

6 μs

xxxxxxxxxx
 
myArray = [2, 6, 5];

4.5 μs

xxxxxxxxxx
 
myDict = Dict("foo" => 3, "bar" => 5, "baz" => 8);

40.9 ms

Arrays are accessed by indexing

xxxxxxxxxx

6 μs

xxxxxxxxxx
 
myArray[1];

2.4 μs

Dictionaries are (unordered) collections of (key, value)-pairs, values are retrieved via corresponding key

xxxxxxxxxx

5.7 μs

xxxxxxxxxx
 
myDict["foo"];

12.7 ms

Loops and control-flow

xxxxxxxxxx

2.9 μs

xxxxxxxxxx
 
my_arr1 = [4, 5, 6];

2.5 μs

xxxxxxxxxx
 
my_arr2 = zeros(length(my_arr1));

7.5 ms

xxxxxxxxxx
 
for i in 1:length(my_arr1)
    my_arr2[i] = 2*my_arr1[i]
end

2.9 ms

xxxxxxxxxx
 
my_arr2;

595 ns

Functions

xxxxxxxxxx

2.9 μs

xxxxxxxxxx
 
function double_vals(my_arr1)
    my_arr2 = zeros(length(my_arr1))
    for i in 1:length(my_arr1)
        my_arr2[i] = 2*my_arr1[i]
    end
    my_arr2
end;

27.8 μs

xxxxxxxxxx
 
double_vals(my_arr2);

11.2 ms

Type inference

xxxxxxxxxx

2.8 μs

Julia is a dynamic language, you don't have to specify types, but you can. Either way, Julia will try to infer types that allow fast subsequent computing.

xxxxxxxxxx

2.6 μs

some_array

Any1

$10$

"foo"

$f a l s e$

xxxxxxxxxx
 
some_array = [10, "foo", false]

25.1 ms

Array{Any,1}

xxxxxxxxxx
 
typeof(some_array)

5.5 μs

xxxxxxxxxx
 
[1, 2, 3];

744 ns

xxxxxxxxxx
 
[true, false, true];

11.4 ms

xxxxxxxxxx
 
[1.0, 2.0, 3.0];

2.1 μs

Abstract and concrete types

xxxxxxxxxx

4.3 μs

Relationships between types

xxxxxxxxxx

10.6 μs

xxxxxxxxxx
 
Int <: Number;

10.1 μs

xxxxxxxxxx
 
Int <: Float64;

6 μs

xxxxxxxxxx
 
Real <: Number;

5.4 μs

xxxxxxxxxx
 
(Array{T, 1} where T) == (Vector{T} where T);

16.3 μs

Parents and children of types

xxxxxxxxxx

7.1 μs

xxxxxxxxxx
 
supertype(Number);

4.5 μs

xxxxxxxxxx
 
subtypes(Number);

117 ms

xxxxxxxxxx
 
subtypes(Real);

52.1 ms

xxxxxxxxxx
 
subtypes(AbstractFloat);

46 ms

Build your own types

xxxxxxxxxx

9.4 μs

xxxxxxxxxx
 
abstract type Person end;

5.1 μs

xxxxxxxxxx
 
struct Professor<:Person
    students::Vector{String}
end

3.1 ms

xxxxxxxxxx
 
Professor <: Person, Person <: Professor;

3.7 μs

xxxxxxxxxx
 
Thomas = Professor(["Matthias", "Maurice"]);

2.4 ms

xxxxxxxxxx
 
typeof(Thomas);

1.3 μs

xxxxxxxxxx
 
Thomas;

591 ns

Multiple Dispatch

xxxxxxxxxx

6.1 μs

Multiple dispatch is the primary paradigm in Julia; it allows to define the same function for multiple type combinations for the input variables. The most appropriate (i.e. fastest) version is then dispatched to a certain input.

xxxxxxxxxx

8 μs

+ is a function in Julia with many different versions (called methods)

xxxxxxxxxx

4.9 μs

xxxxxxxxxx
 
+;

4.2 μs

xxxxxxxxxx
 
+(2, 3) == 2 + 3;

4.1 μs

xxxxxxxxxx
 
@which 1 + 2;

65.4 μs

xxxxxxxxxx
 
@which 1.0 + 2.0;

9.6 ms

xxxxxxxxxx
 
@which 1.0 + 2;

10.7 ms

Create the function mysum with two different methods

xxxxxxxxxx

7.2 μs

xxxxxxxxxx
 
begin
    function my_sum(arr::Vector{Int})
      val = 0
      for i in arr
        val += i
      end
      "sum over ints...", val
    end
​
    function my_sum(arr::Vector{Float64})
      val = 0.0
      for i in arr
        val += i
      end
      "sum over floats...", val
    end
end;

37.2 μs

xxxxxxxxxx
 
my_sum;

506 ns

xxxxxxxxxx
 
my_sum([1, 2, 3, 4]);

7.9 ms

xxxxxxxxxx
 
my_sum([1.0, 2.0, 3.0, 4.0]);

6.2 ms

2. Performance tips

xxxxxxxxxx

3.5 μs

For-loops vs. loop-fusion/vectorisation/broadcasting

Exercise: Which version is the fastest?

xxxxxxxxxx

3.2 μs

xxxxxxxxxx
 
begin
    # x is supposed to be an array like
    x = [i for i in 1:10]
​
    function mult3_v1!(x)
        x .= 3 .* x
    end
​
    function mult3_v2!(x)
        x = 3 * x
    end
​
    function mult3_v3!(x)
        @inbounds for i in 1:length(x)
            x[i] = 3 * x[i]
        end
        x
    end
end;

2.9 ms

xxxxxxxxxx
 
mult3_v1!(x);

47.8 ms

Check box to show solution

xxxxxxxxxx

70.7 μs

You will do it!

xxxxxxxxxx

5.2 ms

In-place functions

Exercise: What could 'in-place' mean? Which version is in-place and what might be faster?

xxxxxxxxxx

3 μs

xxxxxxxxxx
 
begin
    function inplace_v1!(x)
        x .= 3 .* x
        x
    end
​
    function inplace_v2(x)
        y = zeros(Int, length(x))
        y .= 3 .* x
        y
    end
end;

28.1 μs

Check box to show solution

xxxxxxxxxx

60.2 μs

You will do it!

xxxxxxxxxx

8.9 μs

Type-stable code

Exercise: What could 'typle-stable' code mean? Which of the following versions has a type-instability? Which is faster?

xxxxxxxxxx

5.3 μs

xxxxxxxxxx
 
begin
    function divide_v1()
        x = 1
        for i in 1:10
            x = x/2
        end
        x
    end
​
    function divide_v2()
        x = 1.0
        for i in 1:10
            x = x/2
        end
        x
    end
end;

37.9 μs

Check box to show solution

xxxxxxxxxx

59.1 μs

You will do it!

xxxxxxxxxx

8.5 μs

More performance tips:

Use macros and packages to help you find bottlenecks (not guess about runtime, just test it!); such as @btime, @allocated and @profview
Use the @code_warntype macro to find type-instabilities
Copying arrays can be expensive, often Views of arrays are enough and much faster!
Access multi-dimensional arrays along columns since Julia stores in column-major order

More resources:

Official performance tips from the Julia Documentation
A blog post (a bit older, but concepts are still true)

xxxxxxxxxx

4.5 μs

2. Gillespie algorithm

Simple division process, Julia says 'Challenge accepted!'

In the following we have a stochastic simulation algorithm for a simple division process of a single cell with exponential division times (Gillespie algorithm).

Let's check out how Julia code looks like and how fast it is!

xxxxxxxxxx

751 μs

xxxxxxxxxx
 
begin
    λ = 2.0 # division rate
    x0 = 1 # initial cell number
    Δt = 3.0 # time range to simulate
    n = 10000 # number of simulation repeats
​
    function gillespie_alg(x0, λ, Δt, n)
        # preallocate result array
        res = zeros(Int, n)
​
        for i=1:n
            # current cell numbers and current time
            x = x0
            t = 0.0
​
            # simulate
            while true
                # calculate propensity
                prop = x * λ
​
                # draw exponential random time for next event
                τ = rand(Exponential(1.0/prop))
                
                # store single simulation result
                if t + τ > Δt
                    res[i] = x
                    break
                end
                
                # update current cell numbers and time
                x += 1
                t += τ
            end
        end
        # return result
        res
    end
end;

35.2 μs

Let's check if we get expected results! We have a linear Gillespie algorithm, thus we know the analytical result for the mean cell numbers as an ordinary differential equation: $\frac{d x}{d t} = λ \cdot x \Rightarrow x (t) = x_{0} \cdot exp (λ \cdot t)$

xxxxxxxxxx

2.6 μs

xxxxxxxxxx
 
res = gillespie_alg(x0, λ, Δt, n);

80.4 ms

xxxxxxxxxx
 
# expected mean
mean_sol = x0 * exp(λ * Δt);

4.2 ms

xxxxxxxxxx
 
# mean by Gillespie simulations
mean_res = mean(res);

39.1 ms

xxxxxxxxxx

10.6 s

Exercise (optional): Transfer the code to an editor of your choice (Python, R, Matlab) and compare the runtimes with the Julia version!

xxxxxxxxxx

3.1 μs

Check box to show solution

xxxxxxxxxx

70 μs

You will do it!

xxxxxxxxxx

709 μs

3. DifferentialEquations.jl

xxxxxxxxxx

5.1 μs

An ordinary differential equation (ODE) model

Below you find an ODE model for the famous Lotka-Volterra system; an example of pre-ancient theoretical biology (starting around 1910-ish). It describes predator-prey relationships.

xxxxxxxxxx
 
md"""_An ordinary differential equation (ODE) model_
​
Below you find an ODE model for the famous [Lotka-Volterra](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations) system; an example of pre-ancient theoretical biology (starting around 1910-ish). It describes predator-prey relationships.
​
$(PlutoUI.LocalResource("images/pred_prey.png", :width => 650))
"""

1.4 ms

 
function lotka_volterra(du, u, θ, t)
    x, y = u
    α, β, γ, δ = θ
    du[1] = (α - β*y)x # dx =
    du[2] = (δ*x - γ)y # dy =
end;

29.9 μs

α =

β =

γ =

δ =

xxxxxxxxxx

125 μs

xxxxxxxxxx
 
[α, β, γ, δ];

12.1 ms

xxxxxxxxxx
 
begin
    θ = [α, β, γ, δ]
    u0 = [1.0, 1.0]
    prob_lv = ODEProblem(lotka_volterra, u0, (0.0,10.0), θ)
    sol_lv = solve(prob_lv, Tsit5())
end;

1.4 s

xxxxxxxxxx
 
begin
    plot(sol_lv, linewidth=4, xlabel="Time", ylabel="Population size",
        grid=false)
    plot!(size=(350,200))
end

789 ms

Exercise:

Which variable is prey? Which is predator?
What do the parameters α, β, γ and δ mean biologically?

xxxxxxxxxx

4 μs

Check box to show solution

xxxxxxxxxx

70.6 μs

You will do it!

xxxxxxxxxx

9.3 μs

What else does DifferentialEquations.jl offer? Basically everything...

Discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations)
Ordinary differential equations (ODEs)
Split and Partitioned ODEs (Symplectic integrators, IMEX Methods)
Stochastic ordinary differential equations (SODEs or SDEs)
Random differential equations (RODEs or RDEs)
Algebraic differential equations (DAEs)
Delay differential equations (DDEs)
Mixed discrete and continuous equations (Hybrid Equations, Jump Diffusions)
(Stochastic) partial differential equations ((S)PDEs) (with both finite difference and finite element methods)

And a whole ecosystem around it (SciML):

Parameter estimation
Sensitivity analysis
Bayesian inference
Neural networks + Differential equations
...

xxxxxxxxxx

8.5 μs

4. Turing.jl

xxxxxxxxxx

2.8 μs

Bayesian inference with probabilistic programming (Poisson example)

Let's load some data: (two data sets available data5 and data50)

xxxxxxxxxx

1.6 ms

xxxxxxxxxx
 
data_poisson_dict = load("data/data_poisson.jld");

2.1 s

xxxxxxxxxx
 
data_poisson = data_poisson_dict["data50"];

3.5 μs

xxxxxxxxxx

1.7 s

Define a prior:

xxxxxxxxxx

3.2 μs

xxxxxxxxxx
 
begin
    gamma_α = 5
    λ_poi_prior = Gamma(gamma_α, 1) # Uniform(0, 20) # Uniform(0, 5)
end;

7.3 ms

xxxxxxxxxx

637 ms

Set up the Turing model:

xxxxxxxxxx

2.8 μs

xxxxxxxxxx
 
@model function poisson_fit(data_poisson)
    λ ~ λ_poi_prior
    for i in 1:length(data_poisson)
        data_poisson[i] ~ Poisson(λ)
    end
end;

58.9 μs

Run the inference:

xxxxxxxxxx

3 μs

xxxxxxxxxx
 
begin
    iterations = 1000
    ϵ = 0.05
    τ = 10
    chain_poi = sample(poisson_fit(data_poisson), HMC(ϵ, τ), iterations)
end;

11.8 s

xxxxxxxxxx

956 ms

xxxxxxxxxx

1.3 s

xxxxxxxxxx

1 s

Exercise: Change the Gamma prior to an Uniform prior! Which kind of uniform prior is a good choice?

xxxxxxxxxx

2.6 μs

Check box to show solution

xxxxxxxxxx

71.1 μs

You will do it!

xxxxxxxxxx

9.4 μs

Some theoretical background for this example

Using Bayesian inference: How to learn something about a parameter θ of a model M given some data D? We first interpret θ and D as random variables. We then incorporate (subjective) prior knowledge about θ as a probability distribution p(θ) (the 'prior'). Our goal is now to calculate the 'posterior' p(θ|D) (the new information for θ given the data D). To get this we use Bayes' rule: $p (θ | D) = \frac{p (D | θ) \cdot p (θ)}{p (D)}$ (in words, posterior = likelihood x prior / evidence ).

What is probabilistic programming? Probabilistic programming is a programming paradigm in which probabilistic models are specified and inference for these models is performed automatically. It represents an attempt to unify probabilistic modeling and traditional general purpose programming in order to make the former easier and more widely applicable (ref: Wiki).

Constitutive mRNA expression at steady state: The stationary/steady state distribution for mRNA numbers X with linear, constant synthesis and degradation rates γ and δ can be shown to be the Poisson distribution $X \sim Poi (λ)$ with $λ = \frac{γ}{δ}$ . Further, to get an analytical result for the posterior p(λ|D), we will make use of the Gamma distribution as its a conjugate prior to the Poisson likelihood; we have p(λ|D) (Gamma) ∝ p(D|λ) (Poisson) x p(λ) (Gamma).

xxxxxxxxxx

27.3 μs

5. Turing.jl and DifferentialEquations.jl

xxxxxxxxxx

5.5 μs

Simple division process, continued

What we have learned so far:

Principles of Julia, basic syntax, performance tips
How to Gillespie-simulate a simple cell division process
How to run Differential Equations in Julia
How to estimate model parameters using Turing

Now let's bring it all together! Goal of this last part:

Use our Gillespie algorithm to create some in silico data for the division process
Use DifferentialEquations.jl to build an ODE model for the mean cell numbers
Use Turing.jl on top of that to estimate our ODE model parameters, given our Gillespie data!

xxxxxxxxxx

893 μs

1. Create some in silico data (with our 'true' parameter value for λ):

xxxxxxxxxx

7.1 μs

 
begin
    init_cells = 1 # initial cell number
    λrate = 2.0 # division rate
    timespan = 3.0 # time range to simulate
    datapoints = 40 # number of simulation repeats
​
    # Gillespie simulations
    gill_data = gillespie_alg(init_cells, λrate, timespan, datapoints)
end;

351 μs

xxxxxxxxxx

114 ms

2. Build an ODE model:

Exercise: Complete the function `simple_div_ode` (below) with the correct ODE model!

xxxxxxxxxx

16 μs

Check box to show solution

xxxxxxxxxx

75 μs

You will do it!

xxxxxxxxxx

8.2 μs

xxxxxxxxxx
 
function simple_div_ode(du, u, θ, t)
    ### ADD YOUR CODE HERE (replace nothing)
    # nothing
    du[1] = θ[1] * u[1]
end;

19 μs

xxxxxxxxxx
 
begin
    prob_div = ODEProblem(simple_div_ode, [init_cells],
                            (0.0, timespan), [λrate])
    sol_div = solve(prob_div, Tsit5())
end;

2.3 ms

xxxxxxxxxx

58.8 ms

3. Setup the Turing model:

Exercise (for the experts): Complete the Turing model below, assuming a Normal/Gaussian error model for our likelihood!

xxxxxxxxxx

9.1 μs

Check box to show solution

xxxxxxxxxx

73.5 μs

You will do it!

xxxxxxxxxx

9 μs

xxxxxxxxxx
 
# priors 
# (used for plotting, should be identical to priors in our Tuing model)
begin
    λest_prior = Uniform(0.0, 4.0)
    σ_div_prior = Uniform(0, 800)
end;

11.9 ms

xxxxxxxxxx

1.3 s

xxxxxxxxxx

1.1 ms

 
begin
    @model function simple_div_fit(gill_data)
        λest ~ Uniform(0.0, 4.0)
        σdiv ~ Uniform(0, 800)
​
        # save_everystep=false provides only the last time point
        prob = remake(prob_div, p=[λest])
        μdiv = solve(prob,Tsit5(), save_everystep=false)
​
        for i in 1:length(gill_data)
            ### ADD YOUR CODE HERE (replace nothing)
            # nothing
            gill_data[i] ~ Normal(μdiv[2][1], σdiv)
        end
    end
end;

105 μs

4. Bayesian inference:

Enjoy fast Julia code for fitting an ODE model with full Bayesian inference in just under one second! :)

(Note: We use a gradient-based sampler; automatic differentiation is carried through the complete Turing model including our ODE, which makes this fast)

xxxxxxxxxx

15.7 μs

xxxxxxxxxx
 
chain_div = sample(simple_div_fit(gill_data), NUTS(.65), 1000);

12.2 s

xxxxxxxxxx

3.3 ms

xxxxxxxxxx

3.2 ms

xxxxxxxxxx

240 ms

xxxxxxxxxx
 
# plot(chain_div)

440 ns

6. Further Information and References

Some valuable links

Julia Website
Julia Documentation
Julia Discourse Forum (ideal for asking questions!)
A concise Julia tutorial
JuliaPro from JuliaComputing (like a conda installation for Python)
JuliaAcademy (free course material)

Editors/Programming Environments for Julia

Editors

Juno (in Atom) (current standard)
VS Code extension (maybe the future standard)

Notebooks

Jupyter Lab/Notebook (if you have Jupyter already just install IJulia on top)
Pluto.jl (used here)

Packages and Ecosystems

Larger ecosystems

SciML - Open Source Scientific Machine Learning
The Turing Language - Bayesian inference with probabilistic programming
Flux - Pure-Julia approach to machine learning with Julia's native GPU and AD support
JuliaStats - Statistics and Machine Learning made easy in Julia
BioJulia - Bioinformatics and Computational Biology in Julia
JuliaGraphs - Graph modeling and analysis packages for Julia
JuMP - Modeling language for Mathematical Optimization
JuliaIO - Group for a unified IO infrastructure

More particular packages (sometimes part of ecosystems above)

DifferentialEquations - Multi-language suite for high-performance solvers of differential equations
DataFrames - In-memory tabular data in Julia
Distributions - Julia package for probability distributions and associated functions
Zygote - next-gen automatic differentiation (AD) system for source-to-source AD in Julia
KissABC - Pure julia implementation for efficient Approximate Bayesian Computation
NestedSamplers - Implementations of single and multi-ellipsoid nested sampling

Performance packages

BenchmarkTools - A benchmarking framework for the Julia language
ProfileView - Visualization of Julia profiling data

Notebooks/Visualisation

Plots - Powerful convenience for Julia visualizations and data analysis
Pluto - Lightweight reactive notebooks for Julia
IJulia - Julia kernel for Jupyter

xxxxxxxxxx

17.9 μs

xxxxxxxxxx

9.5 μs

Well done, that's all!

xxxxxxxxxx

10.3 μs

Presenter's notes

Time management

Introduction to Julia 18 min.
Performance tips and a Gillespie algorithm 18 min.
Differential Equations in Julia 18 min.
Bayesian Inference and Probabilistic Programming 18 min.
Everything from above combined :) 18 min.
Further Information/References 0 min.
Total = 5 * 18 min. = 90 min.

Example formula $\frac{d x}{d t} = α x$

Example table

Meaning	Variable
Number of people	people
Average number of slices each person eats	avg
Number of slices on a piece of pizza	slices

Example box

xxxxxxxxxx

14 μs

Solution

Yes that is right, that's a lot of pizza! Excellent, you figured out we need to round up the number of pizzas!

xxxxxxxxxx

132 ms

Example input types

a =

b =

c =

d =

e =

f =

xxxxxxxxxx

172 μs

$50$

""

"Click"

$f a l s e$

"one"

"#000000"

xxxxxxxxxx

6.5 μs

Statistical inference with Julia

Group Retreat Heidelberg, Matthias and Maurice, 15. September 2020

Before we start

Content of this Tutorial

1. Introduction: What is Julia?

The Two Language Problem

Basic syntax

Getting help

Some basic data structures

Loops and control-flow

Functions

Type inference

Abstract and concrete types

Multiple Dispatch

2. Performance tips

For-loops vs. loop-fusion/vectorisation/broadcasting

In-place functions

Type-stable code

2. Gillespie algorithm

Simple division process, Julia says 'Challenge accepted!'

Exercise (optional): Transfer the code to an editor of your choice (Python, R, Matlab) and compare the runtimes with the Julia version!

3. DifferentialEquations.jl

Exercise:

4. Turing.jl

Bayesian inference with probabilistic programming (Poisson example)

Exercise: Change the Gamma prior to an Uniform prior! Which kind of uniform prior is a good choice?

Some theoretical background for this example

5. Turing.jl and DifferentialEquations.jl

Simple division process, continued

Exercise: Complete the function simple_div_ode (below) with the correct ODE model!

Exercise (for the experts): Complete the Turing model below, assuming a Normal/Gaussian error model for our likelihood!

6. Further Information and References

Some valuable links

Editors/Programming Environments for Julia

Packages and Ecosystems

Well done, that's all!

Presenter's notes

Exercise: Complete the function `simple_div_ode` (below) with the correct ODE model!